clear
rng default

%% U
load IdSetM_final

% 11: IdSetM_final(:,2)
U_11=[min(IdSetM_final(:,2)); max(IdSetM_final(:,2))];
clearvars -except IdSetM_final U_11

% 12: IdSetM_final(:,3)
U_12=[min(IdSetM_final(:,3)); max(IdSetM_final(:,3))];
clearvars -except IdSetM_final U_11 U_12

% 21: IdSetM_final(:,5)
U_21=[min(IdSetM_final(:,5)); max(IdSetM_final(:,5))];
clearvars -except IdSetM_final U_11 U_12 U_21

% 22: IdSetM_final(:,6)
U_22=[min(IdSetM_final(:,6)); max(IdSetM_final(:,6))];
clearvars -except IdSetM_final U_11 U_12 U_21 U_22

save('U_intervals.mat', 'U_11', 'U_12', 'U_21', 'U_22')


clear
rng default

%% V
load IdSetW_final

% 11: IdSetW_final(:,2)
V_11=[min(IdSetW_final(:,2)); max(IdSetW_final(:,2))];
clearvars -except IdSetW_final V_11

% 12: IdSetW_final(:,5)
V_12=[min(IdSetW_final(:,5)); max(IdSetW_final(:,5))];
clearvars -except IdSetW_final V_11 V_12

% 21: IdSetW_final(:,3)
V_21=[min(IdSetW_final(:,3)); max(IdSetW_final(:,3))];
clearvars -except IdSetW_final V_11 V_12 V_21

% 22: IdSetW_final(:,6)
V_22=[min(IdSetW_final(:,6)); max(IdSetW_final(:,6))];
clearvars -except IdSetW_final V_11 V_12 V_21 V_22

save('V_intervals.mat', 'V_11', 'V_12', 'V_21', 'V_22')


clear
rng default

%% Phi

load IdSetM_final
load IdSetW_final


% 11: IdSetM_final(:,2), IdSetW_final(:,2)
M=unique(IdSetM_final(:,2), 'rows', 'stable'); 
W=unique(IdSetW_final(:,2), 'rows', 'stable');
sm=size(M,1);
IdSet11=cell(sm,1);
for i=1:sm
    IdSet11{i}=M(i)+W(:); %11
end
IdSet11=vertcat(IdSet11{:});
Phi_11=[min(IdSet11); max(IdSet11)];
clearvars -except IdSetM_final IdSetW_final Phi_11

% 12: IdSetM_final(:,3), IdSetW_final(:,5)
M=unique(IdSetM_final(:,3), 'rows', 'stable'); 
W=unique(IdSetW_final(:,5), 'rows', 'stable');
sm=size(M,1);
IdSet12=cell(sm,1);
for i=1:sm
    IdSet12{i}=M(i)+W(:); %12
end
IdSet12=vertcat(IdSet12{:});
Phi_12=[min(IdSet12); max(IdSet12)];
clearvars -except IdSetM_final IdSetW_final Phi_11 Phi_12

% 21: IdSetM_final(:,5), IdSetW_final(:,3)
M=unique(IdSetM_final(:,5), 'rows', 'stable'); 
W=unique(IdSetW_final(:,3), 'rows', 'stable');
sm=size(M,1);
IdSet21=cell(sm,1);
for i=1:sm
    IdSet21{i}=M(i)+W(:); %21
end
IdSet21=vertcat(IdSet21{:});
Phi_21=[min(IdSet21); max(IdSet21)];
clearvars -except IdSetM_final IdSetW_final Phi_11 Phi_12 Phi_21

% 22: IdSetM_final(:,6), IdSetW_final(:,6)
M=unique(IdSetM_final(:,6), 'rows', 'stable'); 
W=unique(IdSetW_final(:,6), 'rows', 'stable');
sm=size(M,1);
IdSet22=cell(sm,1);
for i=1:sm
    IdSet22{i}=M(i)+W(:); %22
end
IdSet22=vertcat(IdSet22{:});
Phi_22=[min(IdSet22); max(IdSet22)];
clearvars -except IdSetM_final IdSetW_final Phi_11 Phi_12 Phi_21 Phi_22

save('Phi_intervals.mat', 'Phi_11', 'Phi_12', 'Phi_21', 'Phi_22')

load U_CS
load V_CS
Phi_11_CS=U_11_CS+V_11_CS;
Phi_12_CS=U_12_CS+V_12_CS;
Phi_21_CS=U_21_CS+V_21_CS;
Phi_22_CS=U_22_CS+V_22_CS;

save('Phi_CS.mat', 'Phi_11_CS', 'Phi_12_CS', 'Phi_21_CS', 'Phi_22_CS')


clear
rng default

%% Core

load IdSetM_final
load IdSetW_final
load U_CS
load V_CS


        M_ii_ik=unique([IdSetM_final(:,2) IdSetM_final(:,3)], 'rows', 'stable'); %U_11 U_12
        M_kk_ki=unique([IdSetM_final(:,6) IdSetM_final(:,5)], 'rows', 'stable'); %U_22 U_21
        [cx1, cx2] = ndgrid(1:size(M_ii_ik,1), 1:size(M_kk_ki,1)); 
        cx1=cx1(:);
        cx2=cx2(:);
        M_ii_kk_ik_ki=[M_ii_ik(cx1,1) M_kk_ki(cx2,1) M_ii_ik(cx1,2) M_kk_ki(cx2,2)]; %U_11 U_22 U_12 U_21
        
        %Women
        W_ii_ki=unique([IdSetW_final(:,2) IdSetW_final(:,3)], 'rows', 'stable'); %V_11 V_21
        W_kk_ik=unique([IdSetW_final(:,6) IdSetW_final(:,5)], 'rows', 'stable'); %V_22 V_12
        [cx1, cx2] = ndgrid(1:size(W_ii_ki,1), 1:size(W_kk_ik,1));
        cx1=cx1(:);
        cx2=cx2(:);
        W_ii_kk_ik_ki=[W_ii_ki(cx1,1) W_kk_ik(cx2,1) W_kk_ik(cx2,2) W_ii_ki(cx1,2)]; %V_11 V_22 V_12 V_21
        
        %Compose
        %A(h,1)+B(:,1)+A(h,2)+B(:,2)-(A(h,3)+B(:,3))-(A(h,4)+B(:,4))
        BB = W_ii_kk_ik_ki(:,1) + W_ii_kk_ik_ki(:,2) - W_ii_kk_ik_ki(:,3) - W_ii_kk_ik_ki(:,4);
        AA = M_ii_kk_ik_ki(:,1) + M_ii_kk_ik_ki(:,2) - M_ii_kk_ik_ki(:,3) - M_ii_kk_ik_ki(:,4); 
        D_interval= [min(BB) max(BB)] + [min(AA) max(AA)];


D_CS=U_11_CS+V_11_CS+...
       U_22_CS+V_22_CS-...
       (U_12_CS+V_12_CS)-...
       (U_21_CS+V_21_CS);


save('D_interval', 'D_interval')
save('D_CS', 'D_CS')


clear
rng default

%% Marital education premia

load IdSetM_final
load IdSetW_final
load U_CS
load V_CS
load PYX_cond
load PXY_cond

% Men
AvU21=PYX_cond(4)*IdSetM_final(:,5) + PYX_cond(5)*IdSetM_final(:,6) - PYX_cond(1)*IdSetM_final(:,2)-PYX_cond(2)*IdSetM_final(:,3);
AvU21_interval=[min(AvU21), max(AvU21)];

% Women
AvV21=PXY_cond(4)*IdSetW_final(:,5) + PXY_cond(5)*IdSetW_final(:,6) - PXY_cond(1)*IdSetW_final(:,2)-PXY_cond(2)*IdSetW_final(:,3);
AvV21_interval=[min(AvV21), max(AvV21)];

% Logit
AvU21_CS=PYX_cond(4)*U_21_CS + PYX_cond(5)*U_22_CS - PYX_cond(1)*U_11_CS-PYX_cond(2)*U_12_CS;
AvV21_CS=PXY_cond(4)*V_12_CS + PXY_cond(5)*V_22_CS - PXY_cond(1)*V_11_CS-PXY_cond(2)*V_21_CS;

save('Av_interval', 'AvU21_interval', 'AvV21_interval')
save('Av_CS', 'AvU21_CS', 'AvV21_CS')

